Path integral Monte Carlo for dissipative many-body systems 
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We address the possibility of performing numerical Monte Carlo simulations for the thermo- 
dynamics of quantum dissipative systems. Dissipation is considered within the Caldeira-Leggett 
formulation, which describes the system in the path-integral formalism through the inclusion of 
an influence action that is bilocal and quadratic in the system's coordinates. At a first sight the 
usual direct approach of discretizing the path integral could seem feasible, but complications arise 
when one tries to introduce a physically meaningful dissipation kernel: in particular its imaginary- 
time dependence turns out to be severely singular and difficult to evaluate analytically, in spite of 
the simple expressions for its Matsubara components. We therefore propose to face the numerical 
problem using Fourier path-integral Monte Carlo, that can be formulated in two different ways: 
transforming the continuous paths and then truncating the high Fourier components (with possible 
improvements upon the truncation procedure), or performing the Fourier transformation upon the 
discretized paths. The latter choice leads to a simpler formulation and allows for a better control of 
the extrapolation to the limit of infinite Trotter number. The method is implemented for a single 
nonlinear particle with Ohmic dissipation and for a (f> 4 chain with Drude-like dissipation. 



I. INTRODUCTION 

In the last decades the interest in quantum dissipation!!] 
has come mainly from the study of mesoscopic systems, 
which have been experimentally fabricated and theoreti- 
cally analyzed. In such systems, the characteristic quan- 
tum effects involve a macroscopic number of particles. 
The sizeable dimension of the devices implies that the 
relevant dynamical variables can couple to a very large 
number of degrees of freedom of the surrounding envi- 
ronment (or dissipation bath): this coupling can be de- 
scribed macroscopically without caring for the details of 
the interaction, and can result in dramatic changes in 
the behavior of the system. For instance, therdissipative 
phase transition in Josephson-junction arraysa (JJA). 

While the classical thermodynamics is unaffected by 
dissipation, its quantum counterpart is substantially 
modified, and it constitutes therefore an ideal field to 
study the genuine interplay between quantum fluctua- 
tions and dissipation, which leads in general to inter- 
esting physics in the regimes of high quantum coupling 
and/or low temperature. 

The issue of evaluating thermodynamic quantities in a 
quantum- dissipative system was recently-fajced by an ex- 
tension of the effective-potential methodEra, that is very 
fruitful in the regime of intermediate quantum coupling. 
However, a more powerful tool is required when the aim 
is to study dramatic effects, as, for instance, the dissipa- 
tive phase transition from superconducting to insulating 
behavior in JJA predicted by mean-field theory. Unfortu- 
nately, a suitable theoretical approach, allowing a faithful 
comparison with the experimental findings in the regime 
of high quantum coupling, is still lacking. 



In this paper, we discuss an efficient path-integral 
Monte Carlo (PIMC) approach can be implemented. 
In Section |l| the basic formalism and the connection 
with the phenomenological description dissipation are re- 
viewed. The cus tomary approach to Monte Carlo is set 
up in Section III, where some difficulties are pointed out; 
this leads us to consider Fourier PIMC IV basically ex- 



tending the standard approach developed by many au- 
thors in the 80ies, involving the transformation to Mat- 
subara components and their truncation by partial aver- 
aging, that by the way leads to a reformulation of the ef- 
fective potential method. We propose a slightly different 
scheme for the numerical computation framework in Sec- 
tion which overcomes some ambiguities of the former. 
Eventually, in Section [Vl] the latter method is applied 
for two reference models: it appears that working with 
Fourier transformed variables, possibly using the knowl- 
edge of the exact quantum harmonic propagator, gives 
reliable results for many-body systems with reasonable 
numerical effort. 



II. PATH-INTEGRAL FOR THE DISSIPATIVE 
SYSTEM 

A. Formalism 

In this paper we consider the study of dissipation ef- 
fects onto the thermodynamics of a quantum system with 
Hamiltonian 



in 



2 



The Caldeira-Leggett (CL)BtJ model considers the sys- 
tem of interest as linearly interacting with a bath of har- 
monic oscillators, whose coordinates can be integrated 
out from the path integral, leaving the CL euclidean ac- 
tion: 



S[q]=^ h ^\™q 2 (u) + V{q(u) 
r pn 



S {nl) [q] 



n 1 2 



+ S^[q] (2) 



4ft jo jo 



du I du' k(u—u') q(u)—q(u') . (3) 



The kernel k(u) depends on the temperature T = 
and is a symmetric and periodic function of the 
imaginary-time u, k{u) — k{~u) — k((3h — u); its func- 
tional form depends on the spectral density of the envi- 
ronmental batha; moreover, it has a vanishing average, 

Jg h duk(u) = 0. Thanks to the last property, one can 
write the nonlocal dissipative action also as 



SW[q] = ^ / du du'k{u~u') q{u)q{v!) . 
2ft Jo Jo 



(4) 



The density matrix elements in the coordinate repre- 
sentation are expressed by Fcynman's path integral as 



p(q",q') = V[q]e 
Jq> 



S[q] 



(5) 



where the path integration is defined as a sum over all 
paths q(u), with u 6 [0,(3h], q(0) = q' and q(@K) = q", 
and the partition function reads 



,S[q] 



(6) 



The usual procedure for the phenomenological identi- 
fication of k(u) consists in comparing its explicit expres- 
sion (in terms of the dynamical variables of the oscillator 
bath) with the analogous expression of the (retarded) 
damping function j(t) one gets in deriving the (classical 
or quantum) Langevin equation of motion from the same 
composite HamiltonianEI, 

mq + m f dt' 7 (t - t') q(t') + V'(q) = f{t) , (7) 



where f(t) is the fluctuating force. The relation that is 
found between k(u) and 7 (t) can be expressed in a sim- 
ple way as a relation between their respective Matsubara 
transform, 

27T7) 

k n = \ du e~ w ™" k(u) , v n = — , (8) 
Jo Pft 



and Laplace transform, 



and reads 



7 ( 2 ) = / dt er zt 7 (i) 
Jo 



k n = K| l{z=\u n \) 



(9) 



(10) 



Here it is apparent that fco = 0, i.e. the 'local part' 
is assumed to be fully included as a quadratic term in 
the potential. The following completeness/orthogonality 
relations have to be taken into account: 



E 



n— — oo 



pn 6(u) , 



(ii) 



where 5(u) = S(u + f3K) is the periodic delta function, 
and its inverse 



PR p iu n ph _ 1 

d u e w "" = : = (3H 5 n0 . (12) 



B. Ohmic and Drude dissipation 

In the above dynamical equation (R) the bath spec- 
tral density is assumed to be such to reproduce the most 
useful phenomenological models, namely: 

i) Ohmic (or Markovian) dissipation. This is charac- 
terized by the absence of memory in the dissipative 
term, and corresponds to assuming a separation of 
time scales: the time scale with which the bath re- 
sponds to changes in the system is much smaller 
than the system's typical times. In this case dissi- 
pation can be described by one constant parameter, 
7: 



7 (t)= 7 ^-0 + ) 



7 (z) = 7 



(13) 



ii) Drude-like dissipation. Here the bath responds on a 
time scale uj~ 1 which is comparable to the system's 
typical times: 



7(*) =7^ D 



7 (z) = 7- 



+ z 



(14) 



Therefore, there are two parameters which describe 
dissipation, the intensity 7 and the response fre- 
quency (or 'spectral width') w D ; for uj d — > 00, i.e., 
fast bath response, the Ohmic form is recovered. 

Note that 7(2) has the dimension of a frequency, while 
k n is a squared frequency. Only the above two cases will 
be considered in what follows; although, of course, the 
actual physics of a problem could give more appropriate 
definitions of r y(t). 



C. Imaginary-time kernel for Ohmic dissipation 

From the above formulas it follows that the relation 
connecting the imaginary-time kernel k(u) with the (as- 
sumed known) Laplace transform 7(2) of the damping 
function 7(4) is 



fc(u) 



1 



OO 

E 



e^ u K| 7 (* = KI) 



(15) 
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The point is that the simple cases above give rather in- 
volute results for k(u). Let us calculate it in the Ohmic 
case (|l3|), invoking a criterion of mean convergence for 
the rcsummation: 



■j *(«) = E w 



n— — oo 



7T / 7TW . 



While it might be useful to note that 



7T / 7T U 



ITU (3k 2 ITU 

o u cot — = — d u In sm — 



so that one has alternative expressions, 

fcN = ^a u cot- = -9 lt lnsm-, 



(17) 



(18) 



one can see that the requirement k = J du k(u) = is 
not satisfied and that to fulfil it one must subtract from 
the expression found - this is the reason why we used the 
tilde in the notation k{u) - the product of a (periodic) 
delta function S(u) by an infinite constant: 



k(u) = k{u) — fcp S(u) 



where 



k 



0h-e 



~. > 27 ire 
du k(u) — cot — — > 



00 



(19) 



(20) 



one can indeed verify that the correct Matsubara trans- 
form \v n \l is obtained thanks to the cancelation of two 
divergences, 



6ti 



du k(u) (e- iv " 



. HVn f m , iru _ 

+ lm L ducot lm e 

27 7re / 7re \ 
— cot — 1 — cos 



* dx 
— cot x e 

7T 



0(e) + i"/is n (-isignrc) = 7 \u n \ 



III. REAL-SPACE PIMC 



(21) 



The standard PIMC approach consists in approximat- 
ing the partition function (^J) by discretizing the paths 
q(u) on a finite mesh. Namely, the imaginary-time in- 
terval [0,f3h] is divided into P slices of finite duration 



e — (3h/P, P being the so called Trotter number. Each 
whole path {q(u), u 6 [0, (3K\} turns into the P discrete 
quantities {qi — q(is)}, with the periodicity condition 
q Q = qp, and the action becomes: 



S F 



P r 

E 

1=1 



mP 



2j3K A 



(qe - qi-if + 75 V(q t ) 



S 



(nl) 



(22) 



m(3 2 Ti 
' 4P 2 



p 

E 

Ll' = \ 



kl-v (q e - qg,) 2 . (23) 



The partition function is approximated by 



( mP 
\2ttH 2 i3 



P/2 



e e 



(24) 



In the standard PIMC procedure the thermodynamic av- 
erages (say, Gp) obtained from this multiple integral are 
evaluated by a stochastic simulation, e.g., the Metropo- 
lis algorithm for configuration sampling; this is to be 
done for large enough values of P, and -the exact result 
G = Goo is estimated by extrapolating^ the calculated 
values Gp. For the discrete kernel kg that approximates 
the singular function fc(tt), it is reasonable to keep a piece- 
wise approximation, namely, for £ ^ 



e (f+|) 

If- 7P 
kt = - Jduk(u) = — 



for large P one has 



cot-(£- 



7T7 



sm ■ 



n£ 
~~P 



k)-cot-{e- 



p 2 



£<p £ 2 



(25) 



(26) 



and for £ = 



fc = - 



1 r e / 2 

e J-e.ll 

1 p/3h-e/2 

£ Je/2 



du k(u) = - 
du k{u 



s/2 f f3h 

-e/2 
2 7 P 7T 

JpW COt 2P' 



du k(u) 



(27) 



ves all 



The choice (|25|) should be preferred to (|26|) since it en- 
sures the exact vanishing of Y2e kg. However, itjs ap- 
parent that ko does not contribute to the action 

The interaction along the Trotter direction invo : 
pairs (which is very bad from the point of view of the 
code efficiency) although it is rapidly decreasing (~ £~ 2 ). 
This suggest the possibility of cutting the interaction be- 
yond, say, R th neighbors (keeping only \£ — £'\ < R); 
a rough calculation can be made assuming that the ki- 
netic term dominates, i.e. that (xg—xg-1) 2 ~ g 2 /(c 2 tP), 
which gives a ratio between the discarded and the in- 
cluded dissipative interaction energy ~ 1/ In R. In any 
case, it turns out that a simulation along these lines re- 
quires to deal with long-ranged summations whose short- 
range part is highly singular; moreover, if one would like 
to consider more physical dissipation kernels, e.g., the 
Drude one, the calculation of k(u) and of kg becomes 
very involute in spite of the simple expression of k n . 
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IV. FOURIER PIMC WITH CONTINUOUS 
IMAGINARY TIME 

In order to overcome the above mentioned difficulties, 
let us try now to face the problem from another point 
of view: since we know as 'initial input' the Matsub- 
ara components of the kernel, k n , it is worth to explore 
the possibility of using the simulation technique based 
on the sampling of Fourier components of the path q(u). 
We will follow the .scheme of Refs. with 
some modificationsEj that seem to improve upon their ap- 
proach when the so called partial averaging is performed. 

The Fourier transform of the closed path q(u), u G 
[0,/m], q(0) = q{j3h), reads: 



q(u) 



E 



n— — oc 



E 



q n {u) 



q n (u) = 2{x n cos u n u + y n sin u n u) 



(28) 



where q n = x n + iy n = q*_ n since q(u) is real, so that 
x n = x_ n and y n = — y_ n . Using the completeness and 
orthogonality relations, Eqs. ( p"l| ) and (|l2|), the inverse 
transform is found to be 



'In 



du 



q(u) e i 



(29) 



and obviously g = g is the average point of the path. 

In terms of the transformed variables the action (^) 
takes the form 

%] = ^ £ (4 + K) k| 2 + r y V(q(u)) , 

n— — oo 

. ( 30) 

which accounts in a simple way for the nonlocal dissipa- 
tive part, at the price of leaving the integral involving the 
potential, whose argument is to be meant as expressed 
as in Eq. (p8|). The path integral (||) for the partition 
function transforms into 



d 2 q n e -^(-»+fe»)l?»l 2 

ox Pi -y ^v( g („))|, (3i) 

where |<j„| 2 = x 2 + ?/ 2 and d 2 q„ = dx n dy n and 
C = 



(32) 



The measure can be easily checked in the free-particle 
nondissipative limit. One can think this expression as 
the Gaussian average of the last exponential: 



Z = C a-** J dq / cxp LJ** ^ V(q(u))\\ , 



with 



(33) 



(34) 



and the nonvanishing moments 
«» = «>= ' 



(35) 



2/?m i/ 2 + k n 
i.e., the n-th component of q(u) has the variance 

(<£(")) = 4 (((^)> cos2 + ((y 2 )) sin 2 !/„«) ee o(,36) 

= 7T 2 l U ■ ( 3? ) 

pm v£ + k n 

A MC simulation based on Eq. ( |33| ) involves a Metropolis 
dynamics for the Fourier coefficients q, x n , and y n , with 
a truncation of the series (pg|), say, at n = P; this should 
correspond to a standard simulation with Trotter number 
P. 

On the other hand, the authors of of Refs. 
[l3] always expand q(u) after subtracting the initial point 
q ee q(0), in a sin-only series, i.e., 



q(u) = g + ^ a ra sin 



irnu 



(38) 



The difference of this choice resembles the one between 
the use of fixed boundary conditions (stationary waves, 
nonuniform amplitude) instead of periodic boundary con- 
ditions (plane waves, uniform amplitude). 



A. Partial averaging 



The partial averagings improves upon the rude trun- 
cation of the Fourier series [for q(u), and basically relies 
upon the Jensen incqualitylla. Look again at Eqs. ( |3~l| ) 
and (|33|): these can be expressed as a superposition of 
uncorrelated Gaussian averages (F({x n , y n }))n upon the 
variables x n and y n , and for anyone of these averages we 
can choose to approximate 

{e F >„ > e< F >« . (39) 

Therefore, choosing to retain (and simulate) the Fourier 
components up to n = P, one can estimate what is left 
over in the exact average; separating the components that 
we want to keep (up to n = P) from those which are to 
be averaged out, i.e., 



with 



q(u) = q p (u) + £ p (u) 



q p (u) = q+^2q n (u) 



(40) 

(41) 
(42) 



£p(«) = E ' 

n=P+l 

one can immediately get 

00 

<£(«)> = E (*(«)) £a ' 

n=P+l 

00 2 00 2 

^ = E E (43) 

n=P+l ' n=P+l n n 
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and apply the Jensen inequality for this part getting the 
approximate (upper bound for the) partition function as 
a Gaussian average ((• ■ -}) p over the finite set of the first 
2P + 1 variables, 



Z = C I dq (( exp 



lit. 



%v r („„)}) 



with an effective potential V P given as the Gaussian 
smearing 



)) Qp on the scale of a p , 



(45) 



where ((£ 2 



a p . What makes this result appeal- 



ing compared to the previous approaches is the fact that 
a p does not depend on u, as it occurs for the 'station- 
ary wave' approach, so one can expect that even in the 
nondissipative case this could be an improvement for 
PIMC coding. Moreover, note that taking the roughest 
approximation, i.e. P = 0, one gets exactly the recipe 
for the effective potential introduced by Feynman: 



C 



-0n 



dq e 



-/3V (q) 



(46) 



where Vo(q) — \V(q + £o))) a is broadened with 



- 2 Y> 1 P n 

a ° 13m ^ v 2 + k n k n ^o 12m 

n—l "■ 



(47) 



while fi — > for k n — * 0. 



B. The variational effective potential 

In view of improving the technique, one can speculate 
whether it is possible to better account for the harmonic 
part, in the spirit of Refs. [li|[l7],|l8]. Let us first review 
how the improved variational approximation arises in the 
present context. The aim is to incorporate a frequency 
term in the Gaussian averages (|35|), i.e., in the variances 
appearing in Eq. (|33|), and since there is an overall inte- 
gration over q, the frequency lu = oj{q) can depend on it. 
Thus we rewrite Eq. ([n]) as follows 



Z = Cjdqfl / d\ n e -M^+^]M 2 



x exp 



1 du 



C I dq (( exp 



tit, 



1 du 



SV(q(u)) 



(48) 



where 



SV(q(u)) = V(q(u)) -jLU 2 (q(u)-q) 2 + /i , (49) 



and now 



2f3m vl + k n +u 2 ' 



(50) 



f 



0h. 



, sinh f , 



with f(q) = /3hu(q)/2; note that the integral over q al- 
ways stays in front of the Gaussian averages, so that any 
quantities 'inside' it can naturally depend on q, and there 
is no need to emphasize this dependence. Then, taking 
the Jensen approximation for all fluctuating components, 
one gets 



Z > C dqe 



-pVe S (q) 



(52) 



(V(q + 0)) - y^ 2 (9") Mq) + Kq) (53) 



V eS (q) = (SV(q(u)) 



ao(q) 



E 



dm z — ' v, 

n— 1 



n + K 

h 

2moj 



coth / - - 



(54) 



Note that the dependence of V e g on u disappears upon 
averaging. We have now to maximize the r.h.s. of 
Eq. (p^), i.e. to minimize the effective potential (|52 
in order to determine us 2 (q). Since 



(55) 



d[i m 

a cancelation occurs and what is left is the known deter- 
mination, 

dV eS l/, /T w//- , 2/^\<9 a 



du; 2 



l((V"(q + m-rnu 2 (q))^ = . (56) 



This concludes the derivation of the effective potential. 
Note that there is no need to introduce the parameter 
w(q) of Refs. 16 , 17|,|l8| and to optimize it. 



C. Improved partial averaging 

In order to retain the exact calculation of the first P 
fluctuation variables, let us split q(u) as in Eqs. ( |4"Tf|42] ) 
and introduce the frequency uj 2 — iu 2 (qQ, ...,q p 



the 



Gaussian variances we want to approximate, i.e., those 
labeled by n = P+l, oo : 



Z = C [dqf[\^ [d 2 q n e 
•> n=l L n J 



-/3m(V„+fc„)|g„ 



>< n 

n=P+l 

x exp 



(3mv 2 



at, 



d q n e 



-/3m i/J+t„+u ! \q n 



du r 



V(q(u)) 
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= C 



, (57) 



where 



5V(q(u))=V(q P (u)+Uu))--u; 2 e p (u)+fi P , (58) 



with 



r P 







Eln ^ + Eln 



v\ + k n + u 



- 1 



P+i 



In order to perform the partial averaging, we take now 
the Jensen approximation for the Gaussian components 
beyond the P-th one, so the relevant variance is 



9 1 



and the approximation reads 

Z>cjdq^ cxp | - ^ V p (q p («)) | ^ , (61) 

with the effective potential 

^(«) = «% + ^)L P -y w2ft f +^ > ( 62 ) 

that actually depends on {qof-,q P } since g = q P {u) is 
In order to optimize uj we must mini- 



(63) 



given by Eq. ( [41 
mize the integral of the effective potential, 



. d p du 

0= d^ Jo T v ^ qAu) ^ 



which, since d[i p jdui 2 — ma p /2, gives 

Jo n 12 



,iu f-^(«p(«) + 6.)U-H^ 



and definitely 



The effective potential can therefore be written as 
V P (q) = (V(q+U} ap - -f {V"(q+U)) ap +V 



= o, 

(64) 
(65) 



(66) 



or, using the differential operator A p = | ct p d 2 , as 

V p (q) = (l~A p )e A - V(q)+fi p (67) 



Eqs. (pOD and (J65J) are self-consistent for any choice of 
the arguments (qo, ...,q p ). 



D. Low-coupling approximation (LCA) 

In the above framework the frequency uj 2 (go 5 • • • ) <7 P ) de- 
pends on all simulated variables and the self-consistent 
Eqs. (60) and (p|) give rise to a considerable complexity, 



even for one degree of freedom; indeed, one should practi- 
cally solve those equations after each MC move except for 
very simple potentials as the quartic one just discussed. 
Some kind of LCA is then necessary; among several pos- 
sibilities, the most reasonable choices for approximating 
(59) uj 2 (q , ...,q p ) with some uj 2 are: 



(i) leaving the only dependence on go = Q by av- 
eraging over the fluctuation coordinates gi,...,g p , 
which leaves the need of tabulating the resulting 
Wo(<jf) = (w 2 (g , g P ))) p ; therefore, choosing to 
insert to also in the first P Gaussian averages in 
order to better describe the resulting probability 
distribution, one has 

= (V"{q + Zo)) = e A V"(q) , (68) 

with A = \otd 2 and the full pure-quantum spread 



9 00 1 



V 2 + kn + 



(69) 



(ii) taking the above value in the minimum, lu 2 = 
uj 2 (q—q m ) (of course, it is also possible to take the 
improved LCA, i.e. the self-consistent HA (SCHA) 
of uj 2 (q)), so that the above self-consistent equa- 
tions are solved only once. 

The first choice reduces the complexity of the self- 
consistent equations to the same one of the approach of 
Refs. [H||l7j , [ii| and can then be used for problems with 
few degrees of freedom, while the latter appears to be nec- 
essary when facing many-body problems. In both cases 
the effective potential has to be expanded in the same 
way. After splitting 



uJ 2 {q , ...,q P ) = uj 2 + Suj 2 , 
where (for simplicity the integral is omitted) 

5w 2 = -e^V"{q P {u))-ul, 
m 

we use dfj, p /du> 2 = ma p /2 in expanding 



(70) 



(71) 



m- 1 e*eV"(q I ,(u))-u;Z 



- Mop + \ P V{q P (u)) - - a Qp UIq , (72) 

where terms of order Sui 4 are neglected, and replacing 
this in the effective potential we get 



V p (q) ~ (l-A p +A 0P )e A -y(g)+ M „ f 



a 0P uj 



(1-6 A p ) e A -+ 5A p V(q) + n np - - a op co 2 , 
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and, neglecting terms of order <5A^, the LCA effective 
potential eventually reads 



V P (q) = e Aop V(q)+n OI 



m 



a OP u 



(73) 



Eventually, the expression for the partition function 
suitable for numerical simulation reads 

Z = C J dq^ exp | - ^ V p (q p («)) | j . (74) 

Other possibilities for a LCA are explored in Ref. 
where the above described approach was also imple- 
mented for the Morse potential. 



V. FOURIER PIMC WITH DISCRETE 
IMAGINARY TIME 

In order to numerically evaluate the integral appearing 
in Eq. ffl), we have seen in Section HI that the standard 



PIMC method divides the imaginary-time interval [0, f3K\ 
into P slices of width e = /3U/P, and that the coordinate 
q(u) turns into the discrete quantities qt — q(le). The 
partition function Z and the other macroscopic thermo- 
dynamic quantities are obtained as the P — > oo extrap- 
olation of Eq. (|24|) and of the estimators generated from 
it. 



As mentioned in Section III, the application of this 



direct PIMC approach to a dissipative system is made 
difficult by the fact that the kernel k(u— u') is explic- 
itly known in terms of its Matsubara transform k n , i.e., 
Eq. (fl0|), rather than in the imaginary-time domaintl. 
In fact, it is given in terms of the Laplace transform 
"f(z) of the damping function j(t) appearing in the phe- 
nomenological Langevin equation (0), and we have seen 
for Ohmic dissipation "f(z) — 7 that this makes k(u— u') 
long-ranged, while for a more realistic Drude dissipation 
k(u—u') becomes very hard to evaluate. 

In the previous Section we realized that the Fourier 
path integral is very convenient as far as the treatment 
of the dissipative nonlocal action is concerned, because 
it enters the relevant expressions trough the (assumed 
known) Matsubara components k n . However, the contin- 
uous imaginary-time approach used there has a general 
drawback (also present in the nondissipative case) aris- 
ing from the appearance of the integral of the potential 
in the last exponent of Eq. ( |3l| ) . 

The alternative we propose here is to start from the 
finite-P expression ( |24| ) of the standard PIMC for the 
partition function and make there a lattice (discrete) 
Fourier transform, changing the integration variables 
from qi to q n by setting: 



p-i 



in/P 



(75) 



so that: 
Z P 



n=l 



P-l 



^Jdqjl[dq 



x exp • 



m(3 



p-l 



E (»P.n + k n) ISnf 



P-l 



£En*+E^ 




(76) 



where C is a temperature- independent normalization and 
k n is as given in Eq. (|l0|). Comparing with the previous 
expression (|3l|), two significant differences appear: firstly, 
the last term (integral of the potential along a path) 
is converted to a well-defined summation that doesn't 
require further approximations; secondly, the kinetic- 
energy term contains the finite-P Matsubara frequencies 



vp,, 



2P . 7m 

m^~p 



(77) 



rather than v n = 2im/j3h, which are approached for 
P— >oo. Thanks to these features the expression we 
got is exactly equivalent to the standard finite-P expres- 
sion (p4|), a property that gives us control onto the ex- 
trapolation of the results to P^oo. 

Estimators for the relevant thermodynamic quantities 
can be obtained in the usual waysOEj; for example, from 
the thermodynamic relation U — —dp hiZ, the following 
estimator for the internal energy is found: 



2mP 2 



~P 



r 2 mkn 



(78) 



For a given potential V(q), it is convenient to devise 
a characteristic energy scale e (e.g., the barrier height 
for a double well potential, the well depth for physical 
potentials that vanish at infinity, etc.) and length scale 
a (such that variations of V comparable to e occur on 
this length scale) and write 



V(q) = ev(q/a) 



(79) 



In this way one better deals with the dimensionless co- 
ordinate x = q/cr. If x m is the absolute minimum of 
v(x), the harmonic approximation (HA) of the system is 
characterized by the frequency luq given by 



ma* 



v"=v"(x m ) ; 



(80) 



the coupling parameter g for the system can be defined as 
the ratio between the HA quantum energy-level splitting 
Hujq and the overall energy scale e, 



9 



h 2 v" 



(81) 



The case of weak (strong) quantum effects occurs when 
g is small (large) compared to 1. It is then easy to make 
use of dimensionless variables only, i.e. to give energies 
in units of e, lengths in units of a, frequencies in units 
of loq, and so on; the reduced temperature is t = l/(e/3), 
the reduced damping intensity is 7 = 7/070- 



8 



We can finally write a dimensionless expression for the 
partition function ( |76| ) (for odd Trotter number P = 
2N + 1): 



Z P = Cfrz / dx \ Y[ da n db. 
J J n=l 
( N r 

E 



x exp ■ 



4v"tP z 



sin 1 K n 

P t 



i ) 



(82) 



where x t = i + 2 £; v =1 [a„ cos ^ + 6„ sin ^p], Jf„ = 
k n /u>Q and we have used the symmetry properties of k(u), 
so that Kp- n = K n . The real Fourier variables x, a n and 
b n are dimensionless; the integrals in Eq. (§2h may be 
numerically evaluated by standard Monte Carlo sampling 
techniques, e.g., the Metropolis one. 



VI. FOURIER PIMC WITH DISCRETE 
IMAGINARY TIME: APPLICATIONS 

A. Single particle in the double-well potential 



As a first application we consider a particle in a quar- 
tic double well potential v(x) = (1 — x 2 ) 2 in presence 
of Ohmic dissipation, i.e., k n = 2Tt(t/g)Tn, where T is 
the damping strength in units of luq] the same model was 
already investigated in Ref. || by mea ns of the effective- 
potential method outlined in Section IV B. In Fig. [l] we 
show the Fourier PIMC results for the average poten- 
tial energy (v(x)) at the strong quantum coupling g = 5, 
for different values of the damping strength. The Monte 
Carlo data reported in the figure represent the extrap- 
olation to P — » oo of the results obtained at P = 17, 
33, 65, and 129. First of all, for the non-dissipative sys- 
tem (r = 0) we observe the perfect agreement between 
the exact results (obtained by numerical solution of the 
Schrodinger equation) and the PIMC data, proving the 
reliability of the PIMC code; for the dissipative model, 
the PIMC data provide a novel reference to chesek the va- 
lidity of the previous effective-potential resultsu. In par- 
ticular, the latter turns out to be reliable at lower and 
lower temperature as the damping strength increases: in- 
deed, this is expected since the coordinate fluctuations 
decrease with T, i.e., the coordinate-dependent quantities 
tend to the classical behavior as an effect of dissipation. 
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FIG. 1: Temperature dependence of the average potential 
energy (v(x)) for the single particle in a quartic double well, 
for g = 5 and different values of the damping strength F. 
Empty symbols are PIMC daAa, lines the predictions from 
the effective potential methods and the filled circles are the 
exact results for V — 0. 



V(q) 



„ M 

— T 

2R ^ 



R 



(84) 



8 = 1 

where v(x) = (1 — x 2 ) 2 /8, Q is the quantum coupling and 
e K and R are the kink energ^-aed-length, respectively, in 
the classical continuum limitllSr lr 1 l. In the above Hamil- 
tonian the number of particles in the chain is M and pe- 
riodic boundary conditions are assumed. The canonical 
variables are such that [qi,pj] — iSij and the harmonic 
excitations of this system have the dispersion relation 

n k = Q s K ^Jl + AR 2 sin 2 |. 

We assume independent baths coupled to each de gree 
of freedom of the chainD, so that for this system Eq. (|82| ) 
is easily generalized as 



A/ N 

Z P = Ct~ Y\_ dXi / JJ da in db 

A 1 « " m 1 



(85) 



where the action reads 

M ( N rmp2 
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B. One-dimensional <f> chain 

Let us now consider a many-body dissipative system, 
namely, the quantum discrete <^> 4 chain, whose Hamilto- 
nian may be written as 



H = e K 



Q 2 R ir-^ „ 2 



V(q) 



(83) 



2RtP 



(86) 



with the coordinates expressed in terms of their Fourier 
components as 
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(87) 
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FIG. 2: (qf } vs temperature for the (j> 4 chain, with Q — 0.2, 
R = 5, J7d = 100 and different values of T. The empty 
symbols are PIMC data (extrapolated for P — > oo) and tlpe 
lines are the predictions from the effective potential methodu. 
r = 0: circles and solid line; F = 20: squares and short- 
dashed line; V — 100: triangles and long-dashed line. The 
inset reports the average (v(qi)). 



and the dimensionless temperature reads t — (/3e K ) _1 . 

The average quantities for the dissipative 4 chain pre- 
sented in the figures have been obtained for periodic 
chains of length (~ 10 2 sites) large enough to be rep- 
resentative of the thermodynamic limit for each set of 
physical parameters and by extrapolating to P — > oo 
the results given by simulations at finite P. A Drudc- 



like spectral density, as introduced in Section II B , was 
assumed for the environmental interaction, so that the 
dissipative kernel reads 



K = — - 



l + Q{lD/(2Trtn) 



(88) 



where the dissipation strength T = j/fl and the cut-off 
frequency f^D = wd are also measured in units of the 
characteristic frequency f2 = Qe K . 

The comparison of our PIMC results with those of 
the effective-potential methodu, shown in Figs. ^ and [|, 
clearly indicates that the predictions of the latter are very 
accurate; this accuracy is preserved for fairly large values 
of the quantum coupling, close to the predicted limits of 
applicability of the effective-potential approximation, as 
it appears in Fig. |] which reports data for Q = 1. 

Moreover, in order to get a reliable thermodynamic 
limit, finite-size effects have to be negligible, i.e., the 
number M of sites must be large enough. In this con- 
dition, reaching high Trotter numbers becomes more and 
more computationally demanding and the extrapolation 
to P — > oo problematic. However, such difficult^-, can 
sometimes be overcome by means of a simple trickE3 de- 
vised to improve the bare Monte Carlo outcomes. Ac- 
cording to Eqs. (38) and (44) of Ref. ^2] any finite-P 
PIMC estimate G(P) of a given thermodynamic quan- 
tity G can be corrected by the known error affecting the 
same quantity for the corresponding SCHA system (of 



FIG. 3: {v(qi)) vs temperature for the 4 chain, with Q = 0.2, 
R = 5, Qd = 100 and different values of V. Symbols and lines 
as in Fig. 0. 
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FIG. 4: (qf) vs temperature for the </> 4 chain, with Q = 1, 
R — 3, = 10 and different values of V. The full symbols are 
PIMC data at finite Trotter number (P — 81 for finite F and 
P = 11 for r — » oo); the empty symbols are the extrapolated 
values for P — > oo. Therliines are the predictions from the 
effective potential methodu. 



course, including the dissipative action). This error, that 
can be calculated in a simple way, is just the difference 
between the 'exact' (P — > oo) SCHA value, an d the 

finite-P SCHA estimate G^l(P). Note that any ther- 
modynamic quantity of interest for a quadratic action in 
presence of dissipation at finite P can be obtained start- 
ing from the density matrix given by Eq. (A14) of Ref. || 
with w = and C and A given by Eqs. (36) and (37) of 
the same reference, with oo replaced by N in the limits 
of the summations. We thus correct the bare PIMC data 
G(P) to the improved values 



G HA (P)=G(P)+ G^-G^(P) 



(89) 



This procedure is shown to be very effective also for dis- 
sipative systems, as shown in Fig. pi where the improved 
estimates for the internal energy (7ha(P) display a very 
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FIG. 5: Internal energy (per site) U vs 1/P 2 for the 4 chain 
with Q = 1, R = 3, Qd — 10, T — 20, at the temperature 
t = 0.2. The full triangles are the bare PIMC results U(P), 
while the empty ones report the harmonically-corrected data 
Uua{P)- The lines are linear fits. 



weak dependence on P, at variance with the bare ones 
U(P), allowing us to correctly extrapolate to P — » oo 
by using much smaller values of P. The relevance of the 
harmonic correction, which fully includes the dissipation, 
increases with the dissipation strength. 

We think that the above formulation of the Fourier 
path-integral Monte Carlo can make it affordable to in- 
vestigate the thermodynamics of quantum many-body 
dissipative systems. The examples we reported testify to 
the power of the method and confirm that the effective- 
potential approachou is valid in the expected parame- 
ter range (weak quantum coupling and/or strong dissipa- 
tion). The further developments involve the implementa- 
tion of the Fourier PIMC procedure beyond the limits of 
the effective-potential method. It is expected that it will 
permit to study the behavior of strongly quantum sys- 
tems in presence of dissipation and thus open the possi- 
bility to approach problems like the dissipative transition 
in Josephson-junction arrays. 
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